# The Effects of Lawn Signs on Vote Outcomes: 
# Results from Four Randomized Field Experiments
# Donald P. Green, Jonathan S. Krasno, Alexander Coppock, Benjamin D. Farrer, Brandon Lenoir, and Joshua N. Zingher

# Pooled Analysis

# Uncomment to set your working directory
# setwd("")

# Uncomment to install required packages
# install.packages(c("stargazer", "xtable", "sandwich", "lmtest", "dplyr", "rmeta", "ggplot2"))

# Load Packages, Data, and Functions --------------------------------------

library(dplyr)
library(rmeta)
library(xtable)
library(ggplot2)
library(stargazer)

load("exp_1.rdata")
load("exp_2.rdata")
load("exp_3.rdata")
load("exp_4.rdata")

source("Lawn_Signs_Source.R")

# Conduct Analysis --------------------------------------------------------

exp_1_fit <- lm(dvote12 ~ condition_factor + dvote10 + dvote08 + dvote06 + pvote08, weights=weights, 
                  data=subset(exp_1, treatable==1))

exp_2_fit <- lm(share ~ condition_factor + registereddems + jenningsvoteshare09 + jenningsvoteshare05, weights=weights, data=subset(exp_2, in_whole_exp==1))

exp_3_fit <- lm(share ~  condition_factor + obama_share_2012 + shannon_share_2009 , weights=weights, data=exp_3)

exp_4_fit <- lm(share ~ condition_factor +
                      GOVP2010 + USPP2008 + GOVP2006 + USPP2004 + GOVP2002 + USPP2000,
                    weights = weights, data=exp_4)

exp_1_fit_r <- commarobust(exp_1_fit)
exp_2_fit_r <- commarobust(exp_2_fit)
exp_3_fit_r <- commarobust(exp_3_fit)
exp_4_fit_r <- commarobust(exp_4_fit)

direct_fx <- c(exp_1_fit_r[2,1], exp_2_fit_r[2,1], exp_3_fit_r[2,1], exp_4_fit_r[2,1])
direct_fx_se <- c(exp_1_fit_r[2,2], exp_2_fit_r[2,2], exp_3_fit_r[2,2], exp_4_fit_r[2,2])
indirect_fx <- c(exp_1_fit_r[3,1], exp_2_fit_r[3,1], exp_3_fit_r[3,1], exp_4_fit_r[3,1])
indirect_fx_se <- c(exp_1_fit_r[3,2], exp_2_fit_r[3,2], exp_3_fit_r[3,2], exp_4_fit_r[3,2])

direct_meta <- meta.summaries(direct_fx, direct_fx_se)
indirect_meta <- meta.summaries(indirect_fx, indirect_fx_se)

results_table <- cbind(sprintf("%.3f", c(direct_fx, direct_meta$summary)),
                       paratheses_appender(c(direct_fx_se, direct_meta$se.summary), 3),
                       sprintf("%.3f", c(indirect_fx, indirect_meta$summary)),
                       paratheses_appender(c(indirect_fx_se, indirect_meta$se.summary), 3))

rownames(results_table) <- c("Experiment 1",
                             "Experiment 2",
                             "Experiment 3",
                             "Experiment 4",
                             "Pooled Results")

##replication table seven
# print.xtable(xtable(results_table), 
#              include.colnames =FALSE, 
#              hline.after = 4, only.contents = TRUE, 
#              file = "meta_results.tex")


# Bayesian Updating Graph -------------------------------------------------

bayes_updater <- function(est_1, se_1, est_2, se_2){
  est_pooled <- weighted.mean(c(est_1, est_2), w = 1/(c(se_1, se_2)^2))
  se_pooled <- sqrt(1/((1/se_1^2) + (1/se_2^2)))
  return(c(est_pooled, se_pooled))
}

lawn_learning <- function(prior_mean, prior_sd){
  update_1 <- bayes_updater(est_1 = prior_mean, se_1 = prior_sd,
                            est_2 = direct_fx[1], se_2 = direct_fx_se[1])
  
  update_2 <- bayes_updater(est_1 = update_1[1], se_1 = update_1[2],
                            est_2 = direct_fx[2], se_2 = direct_fx_se[2])
  
  update_3 <- bayes_updater(est_1 = update_2[1], se_1 = update_2[2],
                            est_2 = direct_fx[3], se_2 = direct_fx_se[3])
  
  update_4 <- bayes_updater(est_1 = update_3[1], se_1 = update_3[2],
                            est_2 = direct_fx[4], se_2 = direct_fx_se[4]) 
  return(c(prior_mean, prior_sd, update_1, update_2, update_3, update_4))
}

scenario_1 <- lawn_learning(prior_mean = 0, prior_sd = 0.05)
scenario_2 <- lawn_learning(prior_mean = 0.05, prior_sd = 0.05)
scenario_3 <- lawn_learning(prior_mean = 0, prior_sd = 0.01)

x<-seq(-0.15,.15,0.0001)

df <- rbind(cbind(x, y = dnorm(x,mean = scenario_1[1], sd = scenario_1[2]), time = "Prior", scenario = "Agnostic"),
            cbind(x, y = dnorm(x,mean = scenario_1[3], sd = scenario_1[4]), time = "Update 1", scenario = "Agnostic"),
            cbind(x, y = dnorm(x,mean = scenario_1[5], sd = scenario_1[6]), time = "Update 2", scenario = "Agnostic"),
            cbind(x, y = dnorm(x,mean = scenario_1[7], sd = scenario_1[8]), time = "Update 3", scenario = "Agnostic"),
            cbind(x, y = dnorm(x,mean = scenario_1[9], sd = scenario_1[10]), time = "Update 4", scenario = "Agnostic"),
            cbind(x, y = dnorm(x,mean = scenario_2[1], sd = scenario_2[2]), time = "Prior", scenario = "Optimist"),
            cbind(x, y = dnorm(x,mean = scenario_2[3], sd = scenario_2[4]), time = "Update 1", scenario = "Optimist"),
            cbind(x, y = dnorm(x,mean = scenario_2[5], sd = scenario_2[6]), time = "Update 2", scenario = "Optimist"),
            cbind(x, y = dnorm(x,mean = scenario_2[7], sd = scenario_2[8]), time = "Update 3", scenario = "Optimist"),
            cbind(x, y = dnorm(x,mean = scenario_2[9], sd = scenario_2[10]), time = "Update 4", scenario = "Optimist"),
            cbind(x, y = dnorm(x,mean = scenario_3[1], sd = scenario_3[2]), time = "Prior", scenario = "Skeptic"),
            cbind(x, y = dnorm(x,mean = scenario_3[3], sd = scenario_3[4]), time = "Update 1", scenario = "Skeptic"),
            cbind(x, y = dnorm(x,mean = scenario_3[5], sd = scenario_3[6]), time = "Update 2", scenario = "Skeptic"),
            cbind(x, y = dnorm(x,mean = scenario_3[7], sd = scenario_3[8]), time = "Update 3", scenario = "Skeptic"),
            cbind(x, y = dnorm(x,mean = scenario_3[9], sd = scenario_3[10]), time = "Update 4", scenario = "Skeptic"))

df_positive <- 
  data.frame(means = c(scenario_1[1], scenario_1[3], scenario_1[5], scenario_1[7], scenario_1[9],
                       scenario_2[1], scenario_2[3], scenario_2[5], scenario_2[7], scenario_2[9],
                       scenario_3[1], scenario_3[3], scenario_3[5], scenario_3[7], scenario_3[9]),
             sds = c(scenario_1[2], scenario_1[4], scenario_1[6], scenario_1[8], scenario_1[10],
                       scenario_2[2], scenario_2[4], scenario_2[6], scenario_2[8], scenario_2[10],
                       scenario_3[2], scenario_3[4], scenario_3[6], scenario_3[8], scenario_3[10]),
             time = rep(c("Prior", "Update 1", "Update 2", "Update 3", "Update 4"), 3),
             scenario = rep(c("Agnostic", "Optimist", "Skeptic"), each = 5))

df_joined <- 
  left_join(data.frame(df), df_positive) %>%
  mutate(prop_positive = pnorm(0, mean = means , sd = sds, lower.tail = FALSE),
         facet = paste0(time, ": N(", sprintf("%.3f", means), ",",sprintf("%.3f", sds), ")",
                        "\n pr(ATE > 0) = ", sprintf("%.3f", prop_positive)),
         x=as.numeric(as.character(x)),
         y=as.numeric(as.character(y)))


agnostic_frame <- filter(df_joined, scenario=="Agnostic")
optimist_frame <- filter(df_joined, scenario=="Optimist")
skeptic_frame <- filter(df_joined, scenario=="Skeptic")


g1 <- 
ggplot(agnostic_frame, aes(x=x, y=y)) + 
  geom_line() + geom_vline(xintercept = 0, linetype="dashed") + 
  facet_grid(scenario~facet) + 
  theme_bw() + 
  geom_ribbon(data=subset(agnostic_frame, x > 0 ),aes(ymax=y),ymin=0,
              fill="lightgray",colour=NA,alpha=0.5) + 
  theme(axis.text.y=element_blank(),axis.ticks.y=element_blank(),
        axis.title.x=element_blank(),
        axis.title.y=element_blank(),legend.position="none",
        strip.background = element_blank()) + 
  coord_cartesian(xlim=c(-.11, .11))

g2 <- 
  ggplot(optimist_frame, aes(x=x, y=y)) + 
  geom_line() + geom_vline(xintercept = 0, linetype="dashed") + 
  facet_grid(scenario~facet) + 
  theme_bw() + 
  geom_ribbon(data=subset(optimist_frame, x > 0 ),aes(ymax=y),ymin=0,
              fill="lightgray",colour=NA,alpha=0.5) + 
  theme(axis.text.y=element_blank(),axis.ticks.y=element_blank(),
        axis.title.x=element_blank(),
        axis.title.y=element_blank(),legend.position="none",
        strip.background = element_blank()) + 
  coord_cartesian(xlim=c(-.11, .11))

g3 <- 
  ggplot(skeptic_frame, aes(x=x, y=y)) + 
  geom_line() + geom_vline(xintercept = 0, linetype="dashed") + 
  facet_grid(scenario~facet) + 
  theme_bw() + 
  geom_ribbon(data=subset(skeptic_frame, x > 0 ),aes(ymax=y),ymin=0,
              fill="lightgray",colour=NA,alpha=0.5) + 
  theme(axis.text.y=element_blank(),axis.ticks.y=element_blank(),
        axis.title.x=element_blank(),
        axis.title.y=element_blank(),legend.position="none",
        strip.background = element_blank()) + 
  coord_cartesian(xlim=c(-.11, .11))


##replication figure 2
# png(file = "Bayesian_Learning_For_Lawn_Signs_updated.png", width = 900, height = 600)
multiplot(g1, g2, g3)
# dev.off()

# Heterogeneous FX --------------------------------------------------------

exp_1 <- 
  within(exp_1,{
    party_share <- scale(pvote08)
  })

exp_2 <- 
  within(exp_2,{
    party_share <- scale(jenningsvoteshare09)
  })

exp_3 <- 
  within(exp_3,{
    party_share <- scale(1-obama_2share_2012)
  })

exp_4 <- 
  within(exp_4,{
    party_share <- scale(1-USPP2008)
  })


exp_1_het <- lm(share ~ condition_factor*party_share + dvote10 + dvote08 + dvote06, weights=weights, 
                data=subset(exp_1, treatable==1))

exp_2_het <- lm(share ~ condition_factor*party_share + registereddems + jenningsvoteshare05, weights=weights, data=subset(exp_2, in_whole_exp==1))

exp_3_het <- lm(share ~  condition_factor*party_share + shannon_share_2009 , weights=weights, data=exp_3)

exp_4_het <- lm(share ~ condition_factor*party_share +
                  GOVP2010 + GOVP2006 + USPP2004 + GOVP2002 + USPP2000,
                weights = weights, data=exp_4)

#replication appendix table d.5

#sink("het_fx.tex")
stargazer(exp_1_het, exp_2_het, exp_3_het, exp_4_het,
          se = makerobustseslist(list(exp_1_het, exp_2_het, exp_3_het, exp_4_het)),
          title="Heterogenous Effects by Past Party Support",
          model.names=F, model.numbers=F, style="apsr",
          column.labels = c("Experiment 1", "Experiment 2", "Experiment 3", "Experiment 4"),
          dep.var.labels= "Vote Share of Advertising Candidate",
          covariate.labels=c("Assigned Lawn Signs", "Adjacent to Lawn Signs", "Past Party Support", "Assigned * Support", "Adjacent * Support", NA),
          omit="dvote10|dvote08|dvote06|registereddems|jenningsvoteshare05|shannon_share_2009|GOVP2010|GOVP2006|USPP2004|GOVP2002|USPP2000", 
          omit.labels = "Covariate Adjustment",
          omit.yes.no=c("yes", "no"),
          omit.stat=c("adj.rsq", "ser", "f"),
          star.cutoffs = c(NA, NA, NA),
          notes= c("Robust standard errors are in parentheses.",
                   "Past Party Support is centered at zero and in standard units.",
                   "Covariates for each experiment are listed in in Tables 3-6."),
          notes.append=FALSE)
# sink()

exp_1_het_r <- commarobust(exp_1_het)
exp_2_het_r <- commarobust(exp_2_het)
exp_3_het_r <- commarobust(exp_3_het)
exp_4_het_r <- commarobust(exp_4_het)

direct_ints <- c(exp_1_het_r[8,1], exp_2_het_r[7,1], exp_3_het_r[6,1], exp_4_het_r[10,1])
direct_ints_se <- c(exp_1_het_r[8,2], exp_2_het_r[7,2], exp_3_het_r[6,2], exp_4_het_r[10,2])

meta_direct <- meta.summaries(d = direct_ints, se = direct_ints_se)

indirect_ints <- c(exp_1_het_r[9,1], exp_2_het_r[8,1], exp_3_het_r[7,1], exp_4_het_r[11,1])
indirect_ints_se <- c(exp_1_het_r[9,2], exp_2_het_r[8,2], exp_3_het_r[7,2], exp_4_het_r[11,2])

meta_indirect <- meta.summaries(d = indirect_ints, se = indirect_ints_se)

##replication heterogeneous effects
##mentioned in text of appendix

meta_direct$summary
meta_direct$se.summary

meta_indirect$summary
meta_indirect$se.summary


